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MATHEMATICAL MODELLING BY SYMBOLIC MATHEMATICAL COMPUTATION - 


A COSMOLOGICAL APPLICATION 
James C. Howard 
Ames Research Center 


SUMMARY 


The advantages of symbolic mathematical computation are most evident in 
problems analogous to those described in this report; that is, the application 
of the method to reconstructing several existing cosmological models and their 
associated trajectory equations. The field equations that govern the trajec- 
tories of bodies in space have, in general, large numbers of terms', with each 
term a complicated mathematical expression. The evaluation of these terms and 
the derivation of the equations of the geodesics that describe the trajec- 
tories of bodies in space, require a substantial amount of algebraic manipu- 
lation and symbolic differentiation. For the models considered, the required 
operations were executed with speed and efficiency on an IBM 360/67 computer. 
For example, in the case of the nonhomogeneous Schwarzschild model, the com- 
puter times required to formulate the field and trajectory equations were 
0.74 and 0.30 minutes, respectively. By mechanization of the procedure in the 
manner described, man hours are saved, the possibility of error is reduced, 
and the scope of the inquiry may be extended. 


INTRODUCTION 


Many problems in mathematical physics involve the formulation of complex 
models of physical systems or processes and the manipulation of large sets of 
nonlinear partial differential equations. The manual performance of these 
operations is time consuming, subject to human error, and in certain cases 
impossible. 

In recent years mathematical manipulation and the formulation of mathemat 
ical models have been facilitated by the use of digital computers equipped 
with formula manipulation compilers (FORMAC) and by the development of com- 
puter programs designed to perform a variety of non-numeric operations 
(ref. 1). In references 2 and 3, it is shown how symbolic mathematical com- 
putation was used to obtain the Christoffel symbols of the first and second 
kinds for 12 orthogonal, curvilinear coordinate systems. In references 4 and 
5 it is shown how an IBM 7094 digital computer, equipped with FORMAC, can be 
used to derive the equations of motion of a particle in any curvilinear coordi 
nate system requested by the user. In reference 6 it is indicated how these 
methods can facilitate the derivation of the Navier-Stokes equations of fluid 
motion and the continuity equation. 


The present report indicates how symbolic mathematical computation can be 
used to formulate a variety of cosmological models. Each model is determined 
by the metric of the Riemannian Space, and the only inputs required are the 
metric coefficients of the fundamental quadratic form. For illustrative pur- 
poses, only spherically symmetric static models are considered. It should be 
emphasized, however, that the method described is equally applicable to less 
restricted models, and nonstatic models can be formulated with equal facility. 
The determination of the geodesics that describe the trajectories of bodies in 
space requires that the appropriate potential functions be known. The rela- 
tivistic analog of Poisson's equation, which in the Newtonian theory connects 
a single gravitational potential function with the density of matter, is a 
relation between the potential functions and the components of the energy 
momentum tensor. In general, this relationship gives rise to 10 nonlinear 
partial differential equations. The solution of these equations then yields 
the potential functions and must precede any attempt to obtain the 
corresponding trajectories. 


ANALYSIS 


Field and Trajectory Equations 

Consider the equation of motion of a particle that is moving under the 
influence of gravitational forces. When the equation is written in the nota- 
tion of the tensor calculus, it assumes the following form (ref. 5): 


m 



+ 



dx-* dx k \ - 

dt dt J a i 


V(j) 


where i,j ,k = 1,2,3 


and 





( 1 ) 


( 2 ) 


In these equations the summation convention is assumed. That is to say, 
if in any term an index occurs twice, the term is to be summed with respect to 
that index for all admissible values of the index. 

In relativistic mechanics, equation (1) is replaced by the following 
trajectory equation (ref. 7): 


d 2 x a ( i \ dx^ dx k \ 
ds 2 + \ jk i ds ds ) 


0 


i,j,k — 1,2, 3,4 


(3) 


2 



( 4 ) 


where the line element ds satisfies the fundamental quadratic form 

ds 2 = g. . dx^ dx“* 

*ij 

The Newtonian theory of gravitation connects a single potential function 
<j> with the density of matter. In this theory, the gravitational potential 
function is required to satisfy Poisson’s equation (ref. 7). 

V 2 <f> = -4 ttp (5) 

At all points of space devoid of matter p = 0 and Poisson’s equation 
reduces to Laplace’s equation 


v 2 (f> = 0 (6) 

The relativistic analog of Poisson's equation is the following tensor 
equation (ref. 8) . 


R. . 
ij 


■=■ Rg. - + Ag. . = -kT. . 
2 6 xj lj 
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By raising indices, the field equations can be written in the alternative 

form 


R 1 . - i 6* R + 6 X A 

■3 2 j j 


-kT 1 


3 


( 8 ) 


Contraction of equation (8) yields 


R - 4A = kT 


(9) 


In regions of space devoid of matter, all the components of the energy 
momentum tensor are zero, and equation (9) simplifies accordingly. In this 
case 

R = 4A (10) 


When this result is substituted in equation (7) , the field equations 
assume the form 


R. . = Ag. . (11) 

13 

Nevertheless, in empty space, the trajectories of bodies moving within 
the solar system correspond with great precision to the simpler field equations 
(ref. 8) . 


R. . 


= 0 


( 12 ) 


where 



( 13 ) 


Rij “(ax^xj log} ^ " 9X “ (ij| + lej) |i«) ' lij) 3x 6 1o S /Tsf) 


j - 112,354 


g = 


’1 3 


(14) 


[ij>k] 1 _ !£ii' 


2 1 Sx 3 


ax 8x 


(15) 


{ij} - S^EU.k] 


(16) 


Equation (12) is the relativistic analog of the Lapace equation. It 
represents 10 nonlinear partial differential equations for the 10 unknown 
functions . Once a set of potential functions g_ satisfying equa- 

tion (8) or (12) is found, the corresponding trajectory equations can be 
formulated. 


COMPUTER APPLICATIONS 


For the purpose of illustrating the modelling capability of symbolic 
mathematical computation, a spherically symmetric static field is assumed. 

This assumption implies that the metric tensors g_ are spherically sym- 
metric and independent of the time. Moreover, the metric tensors must be 
chosen in such a way that the line element will reduce to the special relativ- 
ity form for flat space-time. These considerations led to the adoption of the 
following set of metric tensors for anisotropic space. 


L (x 1 ) 

gn = -e 

g 22 = -C* 1 ) 2 

g 33 = -(x 1 sin x 2 ) 2 


^44 


e M(x!) 


J 


( 17 ) 


4 



where the implicit functions Lfx 1 ) and Mfx 1 ) can be adjusted to account for 
the distortion of space in the presence of matter. The corresponding space- 
time interval is 


ds" 


•[- 


^’W ) 2 


- (x 1 dx 2 )' 


(x 1 sin x z dx d ) 


3 -\ 2 


e M(xl) (dx*) 2 ' 


(18) 


where 


dx 4 = c dt 


and 

x 1 radial coordinate 

x polar angle coordinate 

3 

x azimuthal angle coordinate 

and for convenience the velocity of light is assumed equal to 1. 

If the space is assumed to be isotropic, the metric tensors are modified 
as follows: 


®11 = 
§22 = 
§33 = 
§44 = 

The space-time interval in 


L (x 1 ) 


-e L(x hx 1 ) 2 

-e L ^ 1 \x 1 sin x 2 ) 2 
MCx 1 ) 


J 


the isotropic case is 


(19) 


ds 2 = -e L ^ X ^(dx 1 ) 2 + (x 1 dx 2 ) 2 + (x 1 sin x 2 dx 3 ) 2 J + e M ^ X ^(dx 4 ) 2 (20) 


In order to demonstrate the feasibility of using symbolic mathematical 
computation to obtain different models of the universe, a computer program 1 
was written that required only the postulated metric tensors as inputs. 


lr The computer program and related documentation are available from 
Computer Software Management and Information Center (COSMIC), Barrow Hall, 
University of Georgia, Athens, Georgia, 30601. 
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Applications to Cosmological Models 


Anisotropic model — The field equations and the corresponding trajectory 
equations for this condition can be obtained by using the tensors defined in 
equations (17). With these tensors as inputs to a 360/67 digital computer 
which was programmed to derive models of the universe, the following output 
was obtained. 

The metric coefficients determine the gravitational model being studied. 
In order that each run be identified with the correct inputs, the postulated 
metric coefficients are printed out before the main results. In the case 
under consideration, these have the following values 

G (1 , 1) = -E L * 

G (2 , 2) = -X(l) 2 

G(3 ,3) = -SIN 2 (X(2))X(1) 2 

g( 4j 4) = jP-vm 


The program uses the metric tensor inputs to evaluate the Christoffel 
symbols of the first and second kinds. In order to reduce the amount of out- 
put, the Christoffel symbols of the first kind are not printed out. In terms 
of the system coordinates and the unknown functions L and M, the Christoffel 
symbols of the second kind are 

T(l.l,l) = Cl/2) L C1) * (X(l)) 

T (1 , 1 , 2) = 0 
T(l,l,3) = 0 
T(1 , 1 >4) = 0 
T(l, 2,1) = 0 
T(1 ,2,2) = l/XCl) 

T (1 , 2 ,3) = 0 
T (1 , 2 ,4) = 0 
T(l,3,l) = 0 


6 



T(l, 3, 2) = 0 
T(l,3,3) = 1/X(1) 

TCI, 3, 4) = 0 
T(1 ,4, 1) = 0 
T(l,4,2) = 0 
T(l, 4, 3) = 0 

TCI, 4, 4) = Cl/2)M (i;) - (X(l)) 
TC2,1,1) = 0 
TC2,1,2) = 1/XC1) 

T(2 , 1 , 3) = 0 
TC2.1.4) = 0 

T (2 , 2 , 1) = -E -L ‘ ( ' X( ' 1 ^XC1) 

T(2,2,2) = 0 

TC2,2,3) = 0 

TC2,2,4) = 0 

TC2,3,1) = 0 

T (2 , 3 , 2) = 0 

T (2 , 3 , 3) = COSCXC2))/SINCXC2)) 
T(2, 3 ,4) = 0 
T (2, 4, 1) = 0 
TC2,4,2) = 0 
T C2, 4, 3) = 0 
T (2, 4,4) = 0 
T C3, 1 , 1) = 0 
TC3,1,2) = 0 
T C3, 1 , 3) = 1/XC1) 
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T(3,l,4) = 0 
T(3,2,l} = 0 
T(3,2,2) = 0 

T(3,2,3) = CQS(X(2))/SIN(X(2)) 

T (3,2,4) = 0 

T(3,3,l) = _ e " L ' tXtl ^SIN 2 CX(2))XCl) 
T(3,3,2) = -COS (X (2) ) SIN (X(2) ) 

T(3,3,3) = 0 
T(3,3,4) = 0 
T(3,4,l) = 0 
T (3,4,2) = 0 
T(3,4,3) = 0 
T(3,4,4) = 0 
T(4,l,l) = 0 
T(4,l,2) = 0 
T(4,l,3) = 0 

T(4,l,4) = (1/2)M C1) .(XC1)) 

T (4,2,1) = 0 
T(4,2,2) = 0 
T(4,2,3) = 0 
T(4,2,4) = 0 
T (4,3,1) = 0 
T(4,3,2) = 0 
T (4,3,3) = 0 
T(4,3,4) = 0 

T(4,4, 1) = (1/2)E -L*CX(1)) + M.(X(1)) m (1) 


(X(D) 


T (4,4,2) = 0 
T(4,4, 3) = 0 
T(4,4,4) = 0 


Once the Christoffel symbols of the second kind are known the components 
of the Ricci tensor can be derived. The individual components are 

R(1 , 1) = -L (1) • (X(1))/X(1)-(1/4)M (1) • CX(1))L C1) • (X(l)) 

+(1/4)M (1) . 2 (X(1))+(1/2)M (2) * (X(l)) 

R(1 , 2) = 0 
R(1 ,3) = 0 
R(l,4) = 0 
R(2, 1) = 0 

R(2, 2) = - (1/2)E _L ' • (X(1))X(1)+(1/2)E _L ' (X(1)) M (1) • (X(1))X(1) 

+e -L.(X(1))_ 1 

R(2,3) = 0 
R(2,4) = 0 
R(3, 1) = 0 
R(3, 2) = 0 

R(3,3) = -C1/2)E" L ' C1 ) ^ - (X(1))SIN 2 (X(2))X(1) 

+ (1/2)E _L ' (X(1):) M (1) • (X(1))SIN 2 (X(2))X(1)+E" L ‘ '’siN 2 (X(2) ) 

-SIN 2 (X(2) ) 


R(3,4) = 0 
R(4, 1) = o 
R(4,2) = 0 
R(4,3) = 0 
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R(4,4) = -E L ' CX(1)) * M ' tXtDJMa) . (X(l))/X(l) 

+ ( 1 / 4 ) E - L - ( X a))*M. (X(1)) M (1). (X(1 ,) L (1). CX(1)) 

- ( ./2)E- L - (X(I))tM ' (X(1), M (2) . t XCl)) 

G(I,J) and R(I,J) are both known at this stage of the program; therefore 
the Ricci scalar can be obtained. It is given by the following equation 

R = -2E _L ‘ Vx(1) 2 +2E _L ' ( X 0))lO) * (X(l))/X(l) 

-2E" L * • (X(1))/X(1)-(1/2)E _L * (X(1)) M (1) * 2 (X(1)) 

+ (1/2)E _L * ( X U)) m ( 1) • (X(l) ) L^- 1 ^ • (X(1))-E" L ‘ - (X(l)) + 2/X(l) 2 

The preceding information is next used to obtain the field equations. 

The individual equations are 

ET(1 ,1) = E _L ‘ (X(1)) /X(l) 2 +L (;i:) • (X(1))/(E L ‘ (X<;i):) )X(l) 

-e"L* ( x ( 1 )j l C 1 ). (x(i))/x(i)+e" l ' ( x O)) m CD . (x(i))/x(i) 

-(l/4)M i:i) . 2 (X(l))/(E L ‘ (X(1) b 

+ (1/4)M (1) - (X(l))L (i:) • (X(1))/(E L ' (XC1)) )-C1/2)M (:2) -(X(1))/(E L ' CXC1)) ) 
+ (1/4)E" L ‘ (X(;i)) M C1) . 2 (X(l))-(l/4)E _L ' (X(1)) M (1) • (X(1))L (1 ^ • (X(l)) 

+ (1/2)E _L ' CX(D) m (2) • (X(l))-1/X(l) 2 

ET(1,2) = 0 
ET(1,3) = 0 
ET(1 ,4) = 0 
ET (2,1) = 0 
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ET(2,2) = - (1/2)E L ‘ . (x(i))/x(l)+(l/2)E~ L ‘ ( - X ^ 1 ' ) ^M^ 1 - ) • (X(l))/X(l) 

+ (1/ 4)E _L " (X(1 ^M C1) - 2 (X(1))-(1/4)E _L ‘ VM)*™ . (X(1))L (1) • (X(l)) 

+ (1/2)E _L * • (XU)) 


ET(2,3) = 0 

ET(2,4) = 0 

ET(3, 1) = 0 

ET(3,2) = 0 

ET(3,3) = - (1/2) E _L " (X(1)) L (1) • (X(1))/X(1)+(1/2)E _L * • (X(l))/X(l) 

+ (1/4) E _L ' (X(1)) M (1) • 2 (X(1))-(1/4)E -L ' (X(1 ^M (1) . (X(1))L < - 1) . (X(l) ) 
+(1/2)E" L ‘ (X(1)) M (2) • (X(l)) 

ET(3,4) = 0 

ET(4,1) = 0 

ET(4,2) = 0 

ET (4 , 3) = 0 

ET (4,4) = E _L ' (X(1)) /X(1) 2 -E" L ' CXC 1 ))^ 1 ) • (X(l) )/X (1) -1/X(1) 2 


The choice of an energy momentum tensor completes the specification of 
this type of model. Solution of the resulting equations gives rise to the 
components of the potential function. In the case under consideration, the 
solution yields the unknown functions L(x x ) and MCx 1 ). In terms of the 
postulated metric tensor inputs, the computer derives the equation of the 
trajectories as follows: 

A(l) = -(l/2)L ( ' 1 - ) • (X(1))V(1) 2 +E _L ’ ^ X( - 1 ^- ) X(1)V(2) 2 
+E~ L * < - X ^ 1 - ) - ) SIN 2 (X(2))X(l)V(3) 2 
-(1/2)E _L (X(l)) + M* ( x (i))^(l.) # ( x(i))v(4) 2 
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A(2) = -2V(2)V(l)/X(l)+COS(X(2))SIN(X(2))V(3) 2 
A(3) = -2V(3)V(l)/X(l)-2COS(X(2))V(3)V(2)/SIN(X(2)) 
A(4) = (X(1))V(4)V(1) 


Isotropic model — I£ the universe is isotropic, the line element will 
assume the form of equation (20) * When the corresponding metric tensors 
(eq. (19)) were used as inputs to the computer program, the following output 
was obtained. 

The metric coefficients determine the gravitational model being studied. 
In order that each run be identified with the correct inputs, the postulated 
metric coefficients are printed out before the main results. In the case 
under consideration, these have the following values: 

G (1 , 1) = -E 1 ^™ 

G (2 , 2) = -E L ' < X ™X(1) 2 

G(3,3) = -E L ' tX<:i ^SIN 2 (XC2))X(l) 2 

G (4 ,4) = E M ‘ (X(1)) 


The program uses the metric tensor inputs to evaluate the Christoffel 
symbols of the first and second kinds. In order to reduce the amount of out- 
put, the Christoffel symbols of the first kind are not printed out. In terms 
of the system coordinates and the unknown functions L and M, the Christoffel 
symbols of the second kind are: 

T(l.l.l) = C1/2)L (1) *(X(1)) 

TCI. 1,2) = 0 
T(l, 1,3) = 0 
T(1 , 1 ,4) = 0 
T(1 , 2, 1) = 0 

T(1 , 2 , 2) = 1/X(1)+ (l/2)L (i:i • (X(l)) 

T(l,2,3) = 0 
T(l,2,4) = 0 
T(l,3, 1) = 0 
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T(l,3,2) 

= 

0 

TCI, 3, 3) 

= 

i/xci)+Ci/2)L (1) -CXCD) 

TCI, 3, 4) 

= 

0 

TCI, 4,1) 

= 

0 

TCI, 4, 2) 

= 

0 

TCI, 4, 3) 

= 

0 

TCI, 4, 4) 

= 

C1/2)M (1) -CXC1)) 

TC2,1,1) 

= 

0 

TC2,1,2) 

= 

1/X(1)+C1/2)L ( ' 1 ^ ■ CX ( 1 ) ) 

T(2 , 1 ,3) 

= 

0 

TC2,1,4) 

= 

0 

T(2,2, 1) 

= 

-XCD-C1/2)L (1) . (XCD)XCD 

T C2 , 2, 2) 

= 

0 

TC2,2,3) 

= 

0 

T(2,2,4) 

= 

0 

TC2,3,1) 

= 

0 

T C2 , 3 , 2) 

= 

0 

TC2,3,3) 

= 

C0S(XC2))/SINCXC2)) 

TC2,3,4) 

= 

0 

TC2,4,1) 

= 

0 

TC2,4,2) 

= 

0 

T) 2,4,3) 

= 

0 

TC2,4,4) 

= 

0 

TC3,1,1) 


0 

TC3,1,2) 

= 

0 

TC3,1,3) 

= 

l/X(l) + (l/2)L (i:, .CX(l)) 
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T (3,1,4) = 0 
T(3,2,l) = 0 
T(3,2,2) = 0 

1(3,2, 3) = COS (X(2) )/SIN(X(2)) 

T(3,2,4) = 0 

T(3,3,l) = -SIN 2 (X(2))X(1)-1/2L (1) - (X(1))SIN 2 (X(2))X(1) 2 
T (3 , 3 , 2) = -COS (X(2))SIN(X(2)) 

T(3,3,3) = 0 
T(3,3,4) = 0 
T(3,4, 1) = 0 
T(3,4,2) = 0 
T(3,4,3) = 0 
T(3,4,4) = 0 
T(4,l, 1) = 0 
T(4,l,2) = 0 
T (4 , 1,3) = 0 

T (4, 1 , 4) = (1/2)M (1) - (X(l)) 

T(4 ,2, 1) = 0 
T(4 , 2, 2) = 0 
T(4,2,3) = 0 
T(4,2,4) = 0 
T(4,3,l) = 0 
T(4,3,2) = 0 
T(4,3,3) = 0 
T(4,3,4) = 0 

TC4.4.D = ( 1/2 ) E " L * (X(l))+M. (X( 1 )) M ( 1 ) •(X(l)) 
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T(4,4,2) = 0 
T(4,4,3) = 0 
T (4,4,4) = 0 

Once the Christoffel symbols of the second kind are known, the components 
of the Ricci tensor can be derived. The individual components are: 

R(1 , 1) = L (1) . (X(l))/X(l)-(l/4)M (i:) • (X (1) ) L ( " 1 ^ • (X(1)) + (1/4)M (1) . 2 (X(1)) 

+L (2) - (X(1)) + (1/2)M (2;) - (X(l)) 

R(l,2) = 0 
R(l,3) = 0 
R(l>4) = 0 
R(2,l) = 0 

R(2, 2) = (3/2)L (1) - (X(l))X(l) + (l/2)M (i:) . (X(1))X(1) + (1/4)L (1) - 2 (X(1))X(1) 2 
+ (1/4)M^ 1 ^ • (X(l) ) L • (X(1))X(1) 2 +(1/2)L^ 2 ^ ■ (X(1))X(1) 2 

R(2,3) = 0 
R(2 , 4) = 0 
R(3 , 1) = 0 
R(3,2) = 0 

R(3,3) = (3/2)L (1) - (X(1))SIN 2 (X(2))X(1) + (1/2)M ( " 1 ' ) . (X(1))SIN ,2 (X(2))X(1) 

+ (i/4)l <; i) . 2 (XU)) SIN 2 (X (2)) XU) 2 
+ (1/4) • (X(l))L (;i) - (X(1))SIN 2 (X(2))X(1) 2 

+(1/2) L (2) - (X(1))SIN 2 (X(2))X(1) 2 
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R(3 ,4) = 0 
R(4,l) = o 
R(4,2) = 0 
RC4,3) = 0 

R(4,4) = . E -L-( X ( 1 )) + M.(X(1)) M C1). (X(1) )/ X(1) 
.(l/4)E- L *«W) +M -C X ^)M (li . 2 (X(l)) 

- (1/4)E~ L * . (x(l))L < - 1 - ) • (X Cl) D 

-C1/2)B- L '»( 1 ” + M-(XC 1 ))m( 2 ). (x( i ) ) 

G(I,J) and R(I,J) are both known at this stage of the program; therefore 
the Ricci scalar can be obtained. It is given by the following equation 

R = -4E" L ‘ • (X(1))/X(1)-2E _L ’ • (X(l))/X(l) 

— ( 1/ 2) E - L " (XC1)) L (1) . 2 (X(1)) -(1/2)E _L ‘ (X(i:)) M (i:) . 2 C X ( 1 )) 

- (1/2J E -L ' (X(l))L < - 1 - ) • (X(1))-2E~ L ( x ( 1 )) l (2) # 

_ e -L.(X(1)) m ( 2). (x(1)) 

The preceding information is next used to obtain the field equations. 

The individual equations are 

ET ( 1 , 1 ) = -L (1) - (X(1))/(E L ' (X(1)) )X(1)+2E‘ L ' (X(1)) L (1) • (X(l))/X(l) 

+e -L- ( x (1)) m ( 1). (X(1))/X(1)-C1/4)M (1) . 2 CX(1))/(E L ' (x(1):) ) 

+ (1/4)M^ 1 ^ • (X(1))L (1) • (XC1))/(E L ‘ • CXC1) )/ CE L * (X(1)) ) 

-(1/2)M (2) .(X(1))/(E L ’ (X(1)) )+(1/4)E _L ’ (X(1)) L (1) - 2 CX(1)) 

+ (1/4)E" L ‘ C x ( 1 )) m ( 1 ). 2 (x(1 )) + ( 1 / 4 ) E - L -CX(1 )) m ( 1) . (X(1))L (1 - 1 • (X(l)) 
+E" L * (X(i:)) L (2) • (X(l) )+ (1/ 2)E _L ’ CX(1)) M (2:i • (X(l) ) 
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I 

ET(1,2) = 0 
ET(1 , 3) = 0 
ET(1 ,4) = 0 
ET(2,1) = 0 

ET (2,2) = (1/2)E“ L ' • (X(1))/X(1)+(1/2)E _L ‘ . (X(l))/X(l) 

+ (1/4)E _L " (X( 1 );) M (i:) . 2 (X(l)) + (l/2)E' L ' ( x C 1)) l ( 2 ). ( X (l)) 

+ (1/2)E _L ' (X(1)) M (2) - (X(l)) 

ET (2 , 3) = 0 
ET(2,4) = 0 
ET(3, 1) = 0 
ET (3, 2) = 0 

ET (3 , 3) = (1/2) E' L * ( x (1)) l ( 1) . (X(1))/X(1)+(1/2)E _L ‘ (X(1)) M (1) • (X(l))/X(l) 

+ (1/4)E -L * (X(1)) M <;i) . 2 (X(l)) + (l/2)E" L ' ( x ( 1 )) l ( 2 ) . (X(l)) 

+ (1/2)E" L ‘ • (X(l)) 

ET(3,4) = 0 
ET (4 , 1) = 0 
ET (4 , 2) = 0 
ET (4 , 3) = 0 

ET(4,4) = 2E" L ’ (X(1)) L (i:) - (X(l))/X(l)+ (1/4)E _L ‘ (X(1)) L (1) ? (X(l)) 

+e -L- (X(1)) l (2) . (xq-)-, 

The trajectory equations for the isotropic case are: 

A(l) = -(1/2)L (1) . (X(1))V(1) 2 +X(1)V(2) 2 +(1/2)L (1) - (X(1))X(1) 2 V(2) 2 

+ SIN 2 (X(2))X(l)V(3) 2 +(l/2)L < - 1 - ) • (X (1) )SIN 2 (X (2) ) X(l) 2 V(3) 2 
-( 1/2 )E -1j ’ CX(1))^(1) > (x(!))V(4) 2 
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A(2) = -2V(2)V(1)/X(1)-L (1) - (X(l))V(2)V(l)+COS(X(2))SIN(X(2))V(3) 2 
A(3) = -2V(3)V(1)/X(1)-L (1) - (X(l))V(3)V(l)-2COS(X(2))V(3)V(2)/SIN(X(2)) 
A(4) = -M (1) -(X(1))V(4)V(1) 


Static homogeneous models— In the case of a static homogeneous universe, 
it is evident that coordinates can be chosen so that the line element will 
exhibit spherical symmetry around any desired origin, since all parts of the 
universe are permanently alike. Hence, the line element may be taken in the 
spherically symmetric static form of equation (18). In obtaining this form of 
line element, local irregularities in the gravitational field, which would 
occur in the immediate neighborhood of individual stars or stellar systems, 
are neglected. 

For the system described, it can be shown that the components of the 
energy momentum tensor are: 


ET(1 , 1) = ET(2 , 2) = ET(3,3) = 8Trp Q > 

ET(4,4) = -8irp o 

ET(I , J) =0 for I / J 


( 21 ) 


where p Q and p Q are the pressure and density, respectively, as measured by 
an observer who is at least momentarily at rest with respect to the spatial 
axes. The solution of these equations gives rise to the components of the 
potential function. In the case of the field being considered, the solution 
yields the unknown functions Lfx 1 ) and Mfx 1 ) . 

In order to satisfy the conditions of static homogeneity, it can be shown 
that the implicit functions L(x : ) and Mfx 1 ) are subject to the following 
constraints: If the model is homogeneous, the pressure as measured by a local 

observer will be the same everywhere. Again, owing to the assumed homogeneity 
of the model, the density will be the same everywhere. Moreover, the line 
element must reduce to the special relativity form, for flat space-time, owing 
to the known validity of the special theory in such regions. By imposing 
these conditions, it can be shown that there are only three possibilities for 
a static homogeneous model (ref. 8) : 

M = 0 (22) 

L + M = 0 (23) 

L = M = 0 (24) 

These conditions lead respectively to the Einstein, the de Sitter, and the 

special relativity line elements. 


18 



The Einstein model universe— Substitution from equation (22) in 
equation (18) yields the following metric for a homogeneous model which is not 
isotropic. 

ds 2 = £-e L ^ X ^(dx 1 ) 2 - (x 1 dx 2 ) 2 - (x 1 sin x 2 dx 3 ) 2 + (dx 4 ) 2 J (25) 


If the model were assumed to be homogeneous and isotropic, it would be 
necessary to use equation (20) subject to the constraint equation (22), 

Cosmological considerations led Einstein to consider a universe defined 
by the metric (25). When the metric coefficients were supplied as input to 
the computer program, the following output was obtained. 


The metric coefficients determine the gravitational model being studied. 
In order that each run be identified with the correct inputs, the postulated 
metric coefficients are printed out before the main results. In the case 
under consideration, these have the following values: 

G(l.l) - -E L ‘ (X ™ 

G (2 , 2) = -X(l) 2 

G (3 , 3) = -SIN 2 (X[2))X(1) 2 


G(4,4) = 1 


The program uses the metric tensor inputs to evaluate the Christoffel 
symbols of the first and second kinds. In order to reduce the amount of out- 
put, the Christoffel symbols of the first kind are not printed out. In terms 
of the system coordinates and the unknown functions L and M, the Christoffel 
symbols of the second kind are 


T(1 , 1 , 1) 
T(1 , 1 , 2) 
T(l, 1,3) 
T(l,l,4) 
T(l, 2,1) 
TCI, 2,2) 
TCI, 2, 3) 
TCI, 2, 4) 
T(l, 3,1) 


C1/2)L (1) -CX(1)) 

0 

0 

0 

0 

l/xcn 

0 

0 

0 
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T(l,3,2) = 0 
T(l,3,3) = 1/X(1) 

T(l,3,4) = 0 
T(l,4,l) = 0 
T(l,4,2) = 0 
T(l,4,3) = 0 
T(l,4,4) = 0 
T(2,l,l) = 0 
T(2, 1,2) = 1/X(1) 

T(2,l,3) = 0 
T(2, 1,4) = 0 

T(2,2,l) = -E“ L ' CX(1 ^X(1) 

T (2 , 2 , 2) = 0 
T(2,2,3) = 0 
T(2,2,4) = 0 
T(2,3, 1) = 0 
T (2, 3, 2) = 0 

T(2,3,3) = COS (X(2))/SIN(X(2)) 

T (2,3,4) = 0 

T(2,4,l) = 0 

T (2,4,2) = 0 

T(2,4,3) = 0 

T(2,4,4) = 0 

T (3 , 1 , 1) = 0 

T (3, 1 , 2) = 0 

T(3,l,3) = 1/ X Cl) 
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T(3, 1,4) = 0 


T(3,2,l) = 0 
T(3, 2, 2) = 0 

T(3,2,3) = COS (X (2)) /SIN (X(2) ) 
T(3,2,4) = 0 

T(3,3,l) = -E" L ' (X(i;)) SIN 2 (XC2))X(l) 

T(3 ,3, 2) = -COS (X(2))SIN(X(2)) 

T(3,3,3) = 0 

T(3,3,4) = 0 

T (3 ,4 , 1) = 0 

T(3,4 ,2) = 0 

T (3,4,3) = 0 

T(3,4,4) = 0 

T(4, 1, 1) = 0 

T(4, 1,2) = 0 

T (4 , 1 , 3) = 0 

T (4, 1 ,4) = 0 

T(4, 2, 1) = 0 

T (4 , 2 , 2) = 0 

T(4,2,3) = 0 

T(4,2,4) = 0 

T(4,3,l) = 0 

T(4,3,2) = 0 

T(4,3,3) = 0 

T(4,3,4) = 0 

T(4,4,l) = 0 
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T(4,4,2) = 0 
T (4,4,3) = 0 
T(4,4,4) = 0 


Once the Christoffel symbols of the second kind are known, the components 
of the Ricci tensor can be derived. The individual components are 

R(1 , 1) = -L (1) - (X(l))/X(l) 

R(1 , 2) = 0 
R(1 , 3) = 0 
R(l,4) = 0 
R(2 , 1) = 0 

R(2, 2) = -(l/2)E' L ' (X(1)) L (:i:) .(X(l))X(l) + E' L ' CX(:i):) -l 

R(2,3) = 0 

R(2,4) = 0 

R(3 , 1) = 0 

R(3 , 2) = 0 

R(3,3) = - (1/2)E -L * C x CD) l (1) . (X(1))SIN 2 (X(2))X(1)+E" L ‘ (X(1)) SIN 2 (X(2)) 
-SIN 2 (X (2) ) 

R(3,4) = 0 
R(4,l) = 0 
R(4 , 2) = 0 
R(4 , 3) = 0 
R(4,4) = 0 

G (I , J) and R(I,J) are both known at this stage of the program; therefore 
the Ricci scalar can be obtained. It is given by the following equation 

R = -2E _L ’ (X(1 ^/X(1) 2 +2E" L ‘ * (X(l))/X(l)+2/X(l) 2 
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The preceding information is next used to obtain the field equations. 

The individual equations are 

ET (1 , 1) = e" L ' ^ X ^ 1 ^/X(l) 2 +L ( - 1 - ) • (X(1))/(E L ‘ ( - X< - 1 - ) hx(l) 

_ e -L- (X(1)) l ( 1) . ( X (i))/X(l)-1/X(l) 2 
ET(1 , 2) = 0 
ET(1,3) = 0 
ET (1 , 4) •= 0 
ET(2 , 1) = 0 

ET (2,2) = -(1/2)E' L ‘ (X(;i)) L (1) • CX(1))/X(1) 

ET (2 , 3) = 0 
ET(2,4) = 0 
ET (3 , 1) = 0 
ET(3,2) = 0 

ET (3, 3) = -(1/2)E" L ‘ (X(1)) L (1) • (X(l))/X(l) 

ET (3,4) = 0 
ET (4,1) = 0 
ET (4 , 2) = 0 
ET(4,3) = 0 

ET (4 ,4) = e _L ' (X( ^ 1 ^/X(1) 2 -E" L ‘ • (X(l^)/X(l)-1/X(l) 2 

The equations of the corresponding trajectories are: 

A(l) = E' L ‘ (X(1 ^SIN 2 (X(2))X(l)V(3) 2 -(l/2)L (:i) • (X(l) ) V(l) 2 +E _L ' C x C 1 ) ) x C 1 D V C 2) 2 

A(2) = -2V(2)V(l)/X(l)+COS(X(2))SIN(X(2))V(3) 2 

A(3) = -2V(1)V(3)/X(1) -2C0S (X(2))V(2)V(3)/SIN(X(2)) 

A(4) = 0 
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The de Sitter model — As already indicated, the only other general 
relativistic model that is static and homogeneous is the de Sitter universe. 

In the next section the Schwarzschild model will be considered. It will be 
found to have the same form, although not the same content as the de Sitter 
model. Although the Schwarzschild universe is inhomogeneous, the implicit 
functions Lfx 1 ) and Mfx 1 ) that satisfy its field equations also satisfy 
equation (23). In view of these considerations the de Sitter model will not 
be formulated. 

A nonhomogeneous ease— The Schwarzschild model represents a specially 
important application of relativity theory, since it provides a treatment of 
the gravitational field surrounding the sun. This problem was first studied 
by Schwarzschild in 1916, and the results obtained were used to distinguish 
between the predictions of the Newtonian theory of gravitation and the more 
exact predictions of relativity theory. Since the space surrounding the sun 
is assumed to be devoid of matter, all the components of the energy momentum 
tensor are zero. In this case, the field equations have been shown to satisfy 
equation (12). That is: 


R. . = 0 (26) 

Therefore, the components of the Ricci tensor obtained for the anisotropic 
model and satisfying equation (26) yield the components of the potential 
function for the field surrounding a single attracting mass, which is 
spherically symmetric. 

In terms of conventional mathematical symbolism, the Schwarzschild field 
equations assume the following form: 


(27) 


(28) 


(29) 


(30) 


R. . = 0 for i 4 i 
1 3 

The corresponding trajectory equations are: 


'll 


R 
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■p 


d 2 M 


1 dL 


_dM + 2_/ dM Y 

2 dx^ dx 1 4 dx 1 dx 1 4 \dx* / 


_L dL 

x 1 dx 1 


= 0 


1 ) = 0 


R 33 = sin 2 x 2 I e ^ 


‘M-'tf-i?)]- 

K(;S-3)!- 


1 > =0 


R 44 - e 


L 

1 dL 

dM 

i 

d 2 M 

1 dM 

1 1 

/ \ 2_ l 

dM\ 


4 dx 1 

dx 1 

2 

dx 1 dx 1 

x 1 dx 1 

4 I 

Vdx 1 ) J 


= 0 
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d 2 * 1 B [ 1 dL /dx 1 ^ 2 -L i/dx 2 \ 2 
ds 2 " " 2 dx 1 \ds / \ds ) 


+ 


-L 

e 



_1 M-L dM_ / dx 4 \ 2 

2 dx 1 \ds J 


d 2 x 2 

ds 2 




2 

T 

x 


dx 1 dx 2 
ds ds 


+ sin x 2 



d 2 x 3 

ds 2 

d 2 X 4 

ds 2 


(- 

(- 


_2_ dx]_ dx 3 
x 1 ds ds 


2 cot x 2 


dx 2 dx 3 \ 
ds ds / 


dM dx 1 dx 4 \ 
dx 1 ds ds / 


It is seen that 


(31) 


(32) 


(33) 


(34) 


R 33 = sin 2 x 2 • R 22 

and there are therefore only three equations in L and M. In this connection, 
it should be noted that the 10 equations given by equations (8) or (12) are 
not all independent since, theoretically at least, they would then determine 
completely the metric tensor and would restrict the choice of reference system. 
Therefore, there can be no more than six independent conditions between the 
components of Rjj to permit a free choice of coordinate system in 
four-dimensional space (ref. 8). 

The system of 10 nonlinear partial differential equations 


for the 10 unknown functions g^j is very complicated. The general solution 
of this system is not known. However, for the case considered in this paper, 
it is possible to obtain a closed-form solution. It can easily be deduced 
that 


L = -M 


and 


M 

e = 



-L 

e 


2S 



Hence 


g n = - 1/ [ 1 + (a/x 1 )] 
g 22 = -Cx 1 ) 2 

g 33 = -(x 1 sin x 2 ) 2 
gi +4 = 1 + (a/x 1 ) 

If a = -2m, the metric (35) is consistent with the existence of one 
gravitating mass (m) situated at the origin and surrounded by empty space. 

If the metric tensor inputs (eqs. (17)) consisting of unknown functions 
of x 1 are now replaced by the known functions (35), and the program re-run, 
the trajectory equations are obtained in the following form: 




\ (35) 





(36) 


d 2 x 2 + _2_ dx 1 dx 2 
ds 2 x 1 ds ds 


sin x 2 cos x 2 



0 


(37) 



_2_ dx^_ dx^ + 
x 1 ds ds 


2 cot x 2 


dx 2 dx 3 \ 
ds ds J 


0 


(38) 


d 2 x^ _ a dx 1 dx 4 

ds2 (x 1 ) 2 [l + (a/x 1 )] ds ds 


(39) 


CONCLUSIONS 


Symbolic mathematical computation can facilitate the formulation of 
mathematical models. This has been demonstrated by using the method to recon- 
struct several existing cosmological models and their associated trajectory 
equations. It has been shown that such models can be derived with speed and 
efficiency on present generation computers, provided they are equipped with 
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formula manipulation compilers. For example, in the case of the Einstein and 
de Sitter models, the computer times required to formulate the field and tra- 
jectory equations were 0.66 and 0.32 minutes, respectively. For the nonhomo- 
geneous Schwarzschild model, the corresponding times were 0.74 and 0.30 minutes, 
respectively. In addition to saving man-hours and the errors to which humans 
are prone, the method facilitates the study of a greater variety of models. 
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